SoSe2021
→ Dies kann nur durch eine sorgfältige Planung vorweg erreicht werden, nicht nachdem die Daten erhoben wurden.
Eine ausführliche Erklärung zu diesen Punkten ist in S.E. Lazic (2016): Experimental Design for Laboratory Biologists zu finden.
library(labstats) str(fluoxetine)
'data.frame': 20 obs. of 2 variables: $ dose : int 0 0 0 0 0 80 80 80 80 80 ... $ time.immob: int 182 112 206 170 164 158 165 168 182 97 ...
fluoxetine im labstats Pakettime.immob)aov()mod1 <- aov(time.immob ~ factor(dose), data = fluoxetine) summary(mod1)
Df Sum Sq Mean Sq F value Pr(>F) factor(dose) 3 10420 3473 1.982 0.157 Residuals 16 28043 1753
mod2 <- aov(time.immob ~ dose, data = fluoxetine) summary(mod2)
Df Sum Sq Mean Sq F value Pr(>F) dose 1 10161 10161 6.462 0.0204 * Residuals 18 28303 1572 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lm()mod2b <- lm(time.immob ~ dose, data = fluoxetine) summary(mod2b)
Call:
lm(formula = time.immob ~ dose, data = fluoxetine)
Residuals:
Min 1Q Median 3Q Max
-84.96 -12.50 6.30 19.51 73.04
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 170.44000 14.83681 11.488 1.02e-09 ***
dose -0.25200 0.09913 -2.542 0.0204 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 39.65 on 18 degrees of freedom
Multiple R-squared: 0.2642, Adjusted R-squared: 0.2233
F-statistic: 6.462 on 1 and 18 DF, p-value: 0.02044
anova(lm())anova(mod2b)
Analysis of Variance Table
Response: time.immob
Df Sum Sq Mean Sq F value Pr(>F)
dose 1 10161 10160.6 6.462 0.02044 *
Residuals 18 28303 1572.4
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nachdem das Design des Experiments und die statistische Analyse bestimmt wurden, ergeben sich oft folgende 3 Fragen:
\[Power \propto \frac{ES~\alpha~\sqrt{n}}{\sigma^2}\]
Die Teststärke steigt proportional* mit
→ Durch Umstellung der Gleichung kann jede dieser 5 Größen berechnet werden, wenn die anderen 4 gegeben sind.
*Wie genau die Teststärke sich ändert hängt von der Analyse ab.
power.t.test()power.prop.test()power.anova.test()power.t.test(n = NULL, delta = NULL, sd = 1, sig.level = 0.05,
power = NULL,
type = c("two.sample", "one.sample", "paired"),
alternative = c("two.sided", "one.sided"),
strict = FALSE, tol = .Machine$double.eps^0.25)Um die Stichprobengröße für einen unabhängigen t-Test zu ermitteln, brauchen wir einfach nur die Formel der Teststatistik nach `n umwandeln. Es gilt für gleiche Varianzen und Stichprobengrößen:
\[t=\frac{(\bar{X}_1-\bar{X}_2)}{\sqrt{2\frac{s_p^2}{n}}}=\frac{d}{\sqrt{2\frac{s_p^2}{n}}}~~\Rightarrow~~\sqrt{2\frac{s_p^2}{n}}=\frac{d}{t}~~\Rightarrow~~2\frac{s_p^2}{n}=\frac{d^2}{t^2}~~\Rightarrow~~n = \frac{2t^2s_p^2}{d^2}\]
| Kenngröße | Buchfink | Mönchsgrasmücke |
|---|---|---|
| Mittelwert | 1800km | 3000km |
| Standardabweichung s | ±900km | ±1000km |
| Stichprobengröße n | 20 | 30 |
d <- abs(1800-3000) s_p <- sqrt((900^2 + 1000^2) / 2) (n <- (16*s_p^2) / d^2)
[1] 10.05556
| Kenngröße | Buchfink | Mönchsgrasmücke |
|---|---|---|
| Mittelwert | 1800km | 3000km |
| Standardabweichung s | ±900km | ±1000km |
| Stichprobengröße n | 20 | 30 |
d <- abs(1800-3000) s_p <- sqrt((900^2 + 1000^2) / 2) (n <- (16*s_p^2) / d^2)
[1] 10.05556
power.t.test(power=0.8, sd=s_p, delta=d)
Two-sample t test power calculation
n = 10.9148
delta = 1200
sd = 951.3149
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* groupUnd wie groß müsste die Stichprobe sein, wenn wir bereits einen Unterschied von 600km bzw. 300km statistisch nachweisen wollen?
power.t.test(power = 0.8, sd = s_p, delta = 600)
Two-sample t test power calculation
n = 40.44558
delta = 600
sd = 951.3149
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* grouppower.t.test(power = 0.8, sd = s_p, delta = 300)
Two-sample t test power calculation
n = 158.8158
delta = 300
sd = 951.3149
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* grouppower.t.test(n = 20, sd = s_p, delta = d)
Two-sample t test power calculation
n = 20
delta = 1200
sd = 951.3149
sig.level = 0.05
power = 0.9729651
alternative = two.sided
NOTE: n is number in *each* groupIm vorliegenden Experiment betrug der Stichprobenumfang 20 und 30. Angenommen in beiden Fällen war N = 20, wie hoch ist dann die Teststärke gewesen?
power.t.test(n = 5, sd = s_p, delta = d)
Two-sample t test power calculation
n = 5
delta = 1200
sd = 951.3149
sig.level = 0.05
power = 0.4190796
alternative = two.sided
NOTE: n is number in *each* groupWie hoch wäre die Teststärke bei 5 Vogelmessungen gewesen, um einen Unterschied von 1200km statistisch nachweisen zu können?
Angenommen es soll die Untersuchung des Zugverhaltens wiederholt werden, es können aber nur 8 Vögel jeweils untersucht werden. Welche Differenz könnte dann überhaupt statistisch nachgewiesen werden?
power.t.test(power = 0.8, n = 8, sd = s_p)
Two-sample t test power calculation
n = 8
delta = 1433.285
sd = 951.3149
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group→ Die Differenz müsste mindestens 1433.3 km betragen, um überhaupt einen signifikanten Unterschied zwischen den beiden Arten nachweisen zu können.
# Erstelle df mit versch. Kombinationen df <- expand.grid(n = seq(5, 40, 2), sd = c(500,750,1000, 1250), delta = c(500,1000), power = NA) dim(df)
[1] 144 4
# Powerberechnung für jede der Kombinationen
for (i in 1:nrow(df)) {
df$power[i] <- power.t.test(n = df$n[i], sd = df$sd[i], delta = df$delta[i])$power
}df %>% ggplot(aes(n, power, group = sd, color = factor(sd))) + geom_line() + geom_hline(yintercept = 0.8) + facet_grid(~delta, labeller = label_both)
set.seed(123) # Erzeugung von normalverteilten mittleren Differenzen und SDs delta <- rnorm(5000, 1200, 100) sigma <- rnorm(5000, 951, 100) df <- data.frame(delta = delta, sd = sigma) # Berechnung von N für jede Kombination an ES und SD n <- apply(df, 1, function(x) power.t.test(delta = x[1], sd = x[2], power = 0.8)$n)
Wie hoch sollte N sein, um eine 90%ige Wahrscheinlichkeit zu haben, dass die Power 80% beträgt?
quantile(n, 0.9)
90% 14.91489
→ Unsere ursprüngliche Poweranalyse empfahl eine N von 11.
ppower.anova.test(groups = NULL, n = NULL,
between.var = NULL, within.var = NULL,
sig.level = 0.05, power = NULL)Gewichtsunterschiede von Atlantischem Lachs, der in 4 verschiedenen Typen von Netzkäfigen gezüchtet wurde (N = 24).
df <- salmon %>% group_by(cages) %>%
summarise(mean = mean(weight),
sd = sd(weight))power.anova.test(groups = 4, between.var = var(df$mean), within.var = mean(df$sd)^2, power = 0.8)
Balanced one-way analysis of variance power calculation
groups = 4
n = 2.345194
between.var = 0.7663649
within.var = 0.2335063
sig.level = 0.05
power = 0.8
NOTE: n is number in each groupn <- 2:10 p <- power.anova.test(groups = 4, n = n, between.var = var(df$mean), within.var = mean(df$sd)^2) p
Balanced one-way analysis of variance power calculation
groups = 4
n = 2, 3, 4, 5, 6, 7, 8, 9, 10
between.var = 0.7663649
within.var = 0.2335063
sig.level = 0.05
power = 0.6193382, 0.9544039, 0.9968026, 0.9998332, 0.9999929, 0.9999997, 1.0000000, 1.0000000, 1.0000000
NOTE: n is number in each groupplot(n,p$power)
| Funktion | Poweranalyse für |
|---|---|
pwr.2p.test() |
2 Proportionen (gleiches N) |
pwr.2p2n.test() |
2 Proportionen (ungleiches N) |
pwr.anova.test() |
Balanced one-way ANOVA (gleiches N) |
pwr.chisq.test() |
Chi-Quadrat-Test |
pwr.f2.test() |
Generalisierte Lineare Modelle (GLM) |
pwr.p.test() |
Proportion (1 Stichprobe) |
pwr.r.test() |
Korrelation |
pwr.t.test() |
t-Tests (1-2-Stichproben, unabhängig, gepaart) |
pwr.t2n.test() |
t-Test (2-Stichproben, ungleiches N) |
Bei weiteren Fragen: saskia.otto(at)uni-hamburg.de

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License except for the borrowed and mentioned with proper source: statements.
Image on title and end slide: Section of an infrared satellite image showing the Larsen C ice shelf on the Antarctic Peninsula - USGS/NASA Landsat: A Crack of Light in the Polar Dark, Landsat 8 - TIRS, June 17, 2017 (under CC0 license)